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ABSTRACT 


Rough  surface  scattering  plays  a  erueial  role  in  the  statistics  of  acoustic  propagation 
signals,  especially  at  mid-frequencies  and  higher  (e.g.,  acoustic  communications 
systems).  For  many  years,  the  effeets  of  rough  surfaee  scattering  were  eomputed  using 
simple  models  that  were  applied  in  two  dimensions  (2-D)  only. 

A  prescribed  method  of  computing  2-D  rough  surface  scattering  directly  in  a 
parabolic  equation  model  based  on  the  Split-Step  Fourier  algorithm  was  introduced  by 
Tappert  and  Nghiem-Phu  in  the  mid-1980s.  This  method  has  been  successfully 
implemented  in  various  2-D  parabolic  equation  models,  including  the  Monterey  Miami 
Parabolic  Equation  model.  However,  some  scientifie  researeh  of  more  formal  scattering 
predietions  have  suggested  that  out-of-plane,  three  dimensional  (3-D)  scattering  may  lead 
to  signifieant  disparities  in  the  scattered  field  statisties.  Introdueing  a  hybrid 
implementation  for  the  scattering  effect  in  the  field  transformation  equations  using  a  tri¬ 
diagonal  solution  with  the  Pade  approximant  to  obtain  a  system  of  equations  for 
azimuthal  corrections  will  support  predictions  of  the  effect  of  surface  seattering  on  3-D 
propagation,  which  is  critical  in  evaluating  the  variability  in  underwater  acoustic 
propagation.  Results  of  the  3-D  scattering  calculations  obtained  are  compared  with  the 
output  of  basic  2-D  interface  perturbations  utilizing  the  standard  2-D  approaeh. 
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I.  INTRODUCTION 


In  reality,  all  surfaces  are  rough,  and  it  is  often  important  to  know  how  surface 
roughness  changes  the  properties  of  scattered  acoustic  waves  (SAW)  compared  with 
those  reflected  from  a  plane  surface.  The  extent  to  which  surface  roughness  affects  wave 
scattering  and  finally  the  way  it  propagates  is  therefore  of  great  interest.  Over  the  last  few 
decades  this  problem  has  received  considerable  attention  in  several  branches  of  physics, 
meteorology  and  engineering.  Of  major  interest  has  been  the  study  of  electromagnetic 
wave  scattering  from  rough  surfaces,  which  has  been  crucial  in  designing  many  different 
operational  systems  (e.g.,  radar  systems)  where  the  effects  of  the  roughness  of  the  land 
and  sea  surfaces  on  signals  have  been  studied  (e.g.,  Ogilvy  1987).  Similarly,  considerable 
interest  has  been  concentrated  on  the  interaction  of  sound  waves  with  rough  surfaces 
owing  to  its  significance  in  several  fields.  Rough-surface  effects  on  the  acoustics  field 
have  been  extensively  studied  as  well  (e.g.,  Chimenti  and  Lobmis  1998). 

Acoustics,  generally  speaking,  is  one  of  the  primary  means  by  which  the 
underwater  environment  can  be  sampled,  because  acoustic  waves  are  used  in  the 
underwater  medium  in  much  the  same  way  that  electromagnetic  waves  are  used  in  the 
atmosphere.  The  propagation  of  acoustic  waves  in  the  underwater  environment  is 
important  for  several  areas,  including  navigation,  ocean  observation,  detection  of 
submerged  objects,  and  real-time  underwater  communications.  As  the  ocean  environment 
changes,  the  refraction  and  scattering  of  the  acoustic  field  is  affected.  Over  long  time  and 
space  scales,  the  change  in  the  characteristics  of  the  sea-water  throughout  the  water 
column  causes  the  sound  speed  profile  to  change.  This  in  turn  shifts  the  acoustic  paths 
and  changes  the  overall  capacity  for  the  ocean  to  transmit  acoustic  waves. 

Modulating  a  pressure -release  ocean  surface  over  short  time  and  range  scales  has 
been  found  to  be  more  efficient  in  studying  the  acoustic  paths  in  which  virtually  all 
incident  acoustic  energy  is  reflected  in  the  water  column  (Medwin  and  Clay  1998). 
Consequently,  multi-paths  will  be  generated  as  long  as  the  surface  modulates  over  time 
and  space  (Senne  2012). 
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Certainly,  there  are  numerous  environmental  faetors  that  can  influence  the 
acoustic  wave  propagation  in  the  ocean,  such  as  water  temperature  and  salinity,  as  well  as 
sea  bottom  and  surface  perturbations.  The  impact  of  all  these  factors  on  acoustic  wave 
propagation  are  studied  in  order  to  determine  the  performance  level  of  underwater 
acoustic  systems  and  also  to  develop  techniques  for  using  acoustic  signals  to  measure  the 
physical  parameters  of  the  ocean.  Surface  roughness  is  just  one  among  several 
environmental  parameters  that  can  have  noticeable  and  essential  influence  on  underwater 
acoustic  propagation. 

Biot’s  1957  paper  suggested  one  way  to  approach  the  study  of  the  scattering  of 
acoustic  waves  from  a  rough  surface.  The  idea  was  to  treat  the  surface  as  a  collection  of 
hemispherical  protuberances — or  “bosses.”  He  states: 

A  simple  solution  is  developed  for  the  reflected  waves  on  a  rough  surface 
from  a  simple  harmonic  point  source.  It  is  assumed  that  the  roughness  is 
represented  by  a  distribution  of  hemispherical  bosses  whose  size  and 
mutual  distance  are  small  relative  to  the  wavelength.  It  is  shown  that  under 
these  conditions  the  effect  of  the  roughness  is  equivalent  to  a  boundary 
condition  for  the  wave  equation.  This  boundary  condition  embodies  the 
surface  polarization  and  the  mutual  interaction  of  the  bosses. 

Until  recently,  only  2-D  models  have  been  implemented  in  order  to  compute  and 
illustrate  the  rough  surface  effects  upon  the  acoustic  radials  within  the  field  (e.g., 
Shchegrov  1995).  More  recent  3-D  modeling  of  rough  surface  scattering  has  been  based 
on  ray  models  or  localized  calculations  of  full-field  scattering.  As  full  field 
approximations  provide  more  information  than  ray  path  solutions,  a  parabolic  equation 
(PE)  approximation  is  chosen  to  model  the  acoustic  field  out  to  long  range  with  many 
multipath  contributions  from  bottom  and  surface  scattering. 

In  this  study,  a  three  dimensional  (3-D)  ocean  acoustic  model  has  been  combined 
with  a  2-D  rough  surface  formulation  of  a  parabolic  equation  model  to  simulate  the 
effects  of  sea  surface  roughness  on  acoustic  wave  propagation  with  multiple  surface 
interactions.  Rough  sea  surface  realizations  are  generated  in  the  acoustic  model.  An 
approach  to  achieve  high  resolution  and  accurate  results  while  maintaining  computational 
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efficiency  is  applied.  The  results  of  the  scattering  calculations  are  then  compared  with 
results  from  the  simple  2-D  model. 

Throughout  this  thesis,  the  PE  model  is  used  to  determine  the  acoustic  field  at  a 
variety  of  range  and  depth  locations  in  an  underwater  environment.  The  second  chapter 
reviews  the  past  research  and  the  analytical  tools  developed.  It  also  includes  the 
derivation  of  the  PE  equations  for  the  3-D  scattering  model  approach  and  a  description  of 
its  usage  -  including  the  treatment  of  rough  sea  surfaces.  Concerning  the  numerical 
implementations,  Chapter  III  introduces  the  hybrid  implementation  for  the  scattering  term 
using  a  tri-diagonal  solution  in  addition  to  defining  the  rough  surface  perturbation 
equations  that  are  required  for  the  solution.  Results  from  this  acoustic  model  are 
presented  in  Chapter  IV,  which  includes  an  accounting  of  the  field  study  through 
different  surface  perturbations  by  introducing  live  separate  models  while  analyzing  the 
effects  of  rough  surface  on  acoustic  propagation  for  each  model.  The  PE  model  coupled 
with  each  type  of  3-D  scattering  model  is  compared  against  previous  2-D  model  results  in 
Chapter  IV. 

A,  PROBLEM  STATEMENT 

The  intensity  range  of  the  reflected  waves  can  vary  widely  due  to  the  impact  of 
the  evolving  sea  surface  over  a  short  time  scale  (Ogilvy  1987).  Generally,  acoustic  waves 
scatter  from  surfaces  with  wavelengths  that  are  comparable  to  or  larger  than  the  acoustic 
wavelength.  As  the  acoustic  wave  frequency  increases,  the  associated  acoustic 
wavelength  decreases.  Consequently,  smaller  scales  of  roughness  at  the  boundaries 
become  more  important.  This  is  a  key  point  when  considering  sea  surface  modeling.  Sea 
surfaces  are  made  up  of  a  variety  of  wave  amplitudes  and  wavelengths  that  are 
superimposed.  Depending  on  the  acoustic  frequency  being  used,  several  of  these 
individual  surface  wave  components  may  play  important  roles  in  affecting  the  acoustic 
wave.  Eor  example,  high  frequencies,  upon  which  acoustic  communication  systems 
depend,  are  more  susceptible  to  small  fluctuations  in  the  sea  surface  roughness  (Senne 
2012). 
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The  importance  of  the  3-D  out-of-plane  models  has  been  suggested  by  Karjadi 
(2011),  as  he  states: 

Although  the  present  2-D  acoustic  model  provides  a  simple  and 
computationally  efficient  technique  of  examining  the  fluctuations  of 
arrival  times  and  angles  caused  by  surface  waves,  the  2-D  approach  does 
not  directly  address  the  out-of-plane  scattering  caused  by  surface  waves; 
to  properly  model  these  out-of-plane  scattering  effects,  a  full  3-D  acoustic 
wave  model  is  required. 

Moreover,  Senne  (2012)  has  discussed  the  importance  of  the  realistic  evolving 
rough  sea  surface  by  stating: 

Understanding  acoustic  interactions  in  the  vicinity  of  the  rough  sea  surface 
and  maximizing  the  usability  of  underwater  communication  systems  can 
be  achieved  through  the  modeling  of  acoustic  propagation  within  realistic 
ocean  environments  that  include  time-varying  rough  sea  surfaces. 

B,  PURPOSE  STATEMENT 

The  goal  of  the  present  work  is  to  realize  the  general  effects  of  the  environment 
variability,  such  as  rough  surface  scatter,  on  the  3-D  propagation  which  is  critical  in 
evaluating  the  limitations  of  existing  sonar  prediction  models.  If  the  3-D  effects  are  found 
to  be  significant,  this  suggests  sonar  prediction  models  must  be  adapted  to  account  for 
these  effects. 

In  Chapter  II,  the  selected  2-D  and  3-D  techniques  for  rough  surface  scattering 
calculations  are  introduced  along  with  a  description  of  how  the  two  models  are 
integrated. 

Comparisons  between  this  model  output  and  the  output  of  basic  2-D  interface 
perturbation  models  are  considered  to  show  the  feasibility  of  the  approach  and 
demonstrate  the  ability  of  the  model  to  track  acoustic  wave  propagation  under  fluctuating 
sea  surfaces  from  one  range  step  to  another. 
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II.  THEORY 


A  theoretical  procedure  using  a  full  3-D  PE  model  is  developed  in  an  effort  to 
solve  the  acoustic  propagation  problem  at  the  surface  of  the  ocean  in  the  presence  of 
surface  displacements.  This  theory  is  believed  to  be  an  improvement  over  past  efforts  that 
computed  the  effects  of  rough  surface  scattering  using  simple  models  that  were  applied  in 
2-D  only. 

A,  REVIEW  OF  PAST  RESEARCH  AND  ANALYTICAL  TOOLS 

DEVELOPED 

We  discuss  briefly  here  each  of  the  major  techniques  which  have  been  applied  to 
the  issue  of  rough  surface  scattering.  Due  to  the  subject  importance  to  the  various 
scientific  fields,  there  are  numerous  papers  related  to  rough  surface  scattering.  We  have 
selected  a  representative  sample  in  order  to  introduce  previous  techniques  used  and  the 
analytical  tools  developed. 

Probably  the  first  reported  work  on  scattering  from  rough  surfaces  was  by 
Rayleigh  (1878)  who  discussed  the  scattering  of  a  normally  incident  sound  wave  by  a 
corrugated  (sinusoidal)  surface  separating  two  acoustic  media.  He  assumed  that  the 
scattered  field  could  be  written  as  a  sum  of  plane  waves  travelling  away  from  the  rough 
surface.  This  approach  has  since  been  extended  by  other  authors  to  non-normal  incidence 
and  to  random  rough  surfaces  (Ogilvy  1987). 

In  the  perturbation  approach,  the  rough  surface  is  regarded  as  a  perturbation  to  a 
smooth  plane,  and  the  consequent  change  in  the  scattering  coefficient  which  measures 
how  much  of  the  incident  acoustic  waves  are  scattered,  due  to  the  presence  of  roughness, 
is  calculated.  This  approach  requires  that  the  height  deviation  of  the  surface  away  from 
the  smooth  plane  is  everywhere  small  compared  to  the  wavelength  of  the  incident  wave. 
In  addition,  the  slope  of  the  surface  must  be  small  in  comparison  to  unity  (Gilbert  and 
Knopoff  1960). 
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Undoubtedly  the  approaeh  most  prevalent  in  the  literature  is  the  tangent  plane,  or 
Kirehhoff,  method.  Probably  the  first  use  of  this  theory  was  by  Eekart  (1953),  applied  to 
aeoustie  wave  scattering,  and  Davies  (1955)  for  electromagnetic  wave  scattering.  Here 
the  scattering  surface  is  assumed  to  be  everywhere  sufficiently  smooth  so  that  the 
reflection  properties  of  any  point  on  the  surface  are  those  of  a  plane  passing  through  that 
point,  parallel  to  the  local  tangent  to  the  surface.  This  method  leads  to  a  restriction  on  the 
rate  of  change  of  the  surface  gradient  but  no  direct  restriction  on  the  magnitude  of  the 
surface  height  or  surface  gradient  (Ogilvy  1987). 

A  number  of  studies  have  considered  the  mechanisms  governing  interactions 
between  underwater  acoustics  and  rough  sea  surfaces.  When  acoustic  waves  interact  with 
a  perfectly  flat  sea  surface,  they  reflect  from  that  surface  just  as  light  would  reflect  from  a 
mirror.  That  is  to  say,  all  reflected  energy  is  directed  in  the  specular  direction.  As  sea 
surfaces  roughen,  the  acoustic  waves  interact  with  the  varying  scales  of  surface 
roughness,  and  as  a  result  scatter  in  a  multitude  of  directions.  Acoustic  energy  is  thereby 
spread  into  non-specular  reflections.  The  exact  reflection  parameters  are  dependent  on  the 
relationship  between  the  length  scales  of  surface  roughness  and  the  wavelengths  of  the 
impacting  acoustic  waves  (Ogilvy  1987). 

The  exact  solution  of  acoustic  waves  scattering  from  a  rough  sea  surface  can  be 
determined  if  the  incident  acoustic  wave  is  considered  to  be  a  plane  wave,  while  the 
surface  is  assumed  to  be  infinitely  periodic  (Holfford  1981). 

Real  ocean  environments  consist  of  a  multitude  of  surface  wavelengths  and  are 

not  typical  of  these  idealized  environments.  Therefore,  a  variety  of  scattering  models 

have  been  utilized  in  the  literature  to  approximate  the  interactions  between  acoustic 

waves  and  rough  sea  surfaces.  Historically  speaking,  the  Kirehhoff  and  perturbation 

methods  have  been  the  most  prevalent.  The  Kirehhoff  method  is  used  whenever  the 

acoustic  wavelength  is  smaller  than  the  radius  of  curvature  of  the  surface.  This  limits  the 

Kirehhoff  approximation  to  being  primarily  applied  to  high  acoustic  frequencies.  Because 

the  Kirehhoff  formulation  treats  all  points  of  the  scattered  field  on  the  rough  surface  as 

being  locally  flat,  the  approximation  provides  a  series  of  reflections  from  different  points 

along  the  surface.  The  perturbation  method  assumes  the  surface  is  only  slightly  rough, 
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and  treats  the  rough  surfaee  as  a  group  of  perturbations  from  a  smooth  surfaee.  As  sueh,  it 
is  limited  by  the  allowable  magnitudes  of  surface  displacement  and  surface  slope. 
Whereas  the  Kirchhoff  method  is  used  primarily  for  high  acoustic  frequencies,  the 
perturbation  method  is  used  mostly  for  low  acoustic  frequencies,  whose  wavelengths  are 
generally  large  in  comparison  with  the  surface  height  (Ogilvy  1987). 

More  recently,  the  small-slope  approximation  has  gained  in  prominence.  The 
small  slope  approximation  was  analyzed  by  Thorsos  and  Broschat  (1995)  as  a  series  in 
generalized  surface  slope  for  the  pressure-release  (Dirichlet)  boundary  condition.  This 
method  was  later  shown,  in  its  lowest  order,  to  reduce  to  the  perturbation  or  Kirchhoff 
approximations  within  each  of  their  respective  regimes  (Thorsos  and  Broschat  1997). 

To  predict  the  impulse  responses  of  high  frequency  acoustics,  rough  surface 
scattering  is  most  commonly  treated  through  the  application  of  the  previously  mentioned 
scattering  models  to  acoustic  ray  codes.  The  main  benefit  of  these  ray  codes  is  that  they 
provide  a  faster  calculation  scheme  than  full-field  models,  such  as  PE  solutions.  A  lot  of 
work  has  been  done  in  order  to  examine  the  effectiveness  of  a  few  surface  scattering 
mechanisms  applied  to  ray  approximation  models  in  predicting  acoustic  returns  and 
found  that  while  the  small  slope  approximation  gave  the  best  results,  several  of  the 
methods  yielded  inaccurate  results  when  surface  roughness  increased  (Williams,  Thorsos, 
and  Elam  2004). 

PE  models,  as  an  alternative  to  ray-based  models,  have  long  been  used  by  the 
underwater  acoustics  community  to  study  sound  propagation.  By  maintaining  full-field 
characteristics,  PE  models  are  able  to  track  both  the  real  and  imaginary  amplitude  and 
phase  values  of  the  wave  field  over  the  entire  water  column,  giving  a  total  field  at  all 
depths  and  ranges.  PE  models  are  also  able  to  track  fluctuations  in  the  acoustic  channel 
over  short  time  steps,  which  enable  the  calculation  of  coherence  between  successive 
model  runs.  Implementing  rough  surfaces  into  PE  models  follows  a  different  approach 
than  the  ray  codes.  Rather  than  calculate  the  scattering  statistics  of  the  rough  surface  and 
then  use  that  information  to  alter  the  acoustic  wave,  PE  models  instead  adapt  the  range¬ 
marching  scheme  to  incorporate  rough  surfaces  into  the  acoustic  field  (Senne  2012). 
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The  use  of  PE  models  for  acoustic  frequencies  that  are  typical  of  communications 
signals  has  historically  been  sparse  due  to  the  computational  effort  required.  However, 
recent  advances  in  computer  systems  have  allowed  for  faster  calculations  that  enable  PE 
methods  to  be  a  viable  option  even  at  moderate  to  high  acoustic  frequencies  (Senne 
2012). 

B,  DERIVATION  OF  PE  EQUATIONS  FOR  THE  3-D  SCATTERING 
MODEL 

Eollowing  the  development  of  Smith  (2012),  we  begin  with  a  3-D  version  of  the 
standard  PE  approximation  in  Cartesian  coordinates, 

=  .  ^>ri(x.y)  (1) 

where  'P  is  the  PE  “field  function,”  and  T  and  U  are  the  “kinetic”  and  “potential” 

energy  operators  of  the  PE  Hamiltonian.  Using  the  Standard  Parabolic  Equation  (SPE) 
approximation  (e.g..  Smith  2001),  these  operators  are  defined  in  three  dimensions  as 

1  ^2  1  o2  1 

Tov  =  T.dspe  = - - .  and  U  =  =  — («'  - 1).  (2) 

op  iDSPE  2kl  dy^  2 

The  rough  surface  boundary  condition  is  introduced  by  defining  the  location  of 
the  pressure  release  rough  surface  at  depth 

z  =  7(x,y).  (3) 

The  pressure  release  condition  then  requires 

'P(z  =  7(x,y),x,y)  =  0.  (4) 

Introducing  an  image  field  above  the  surface  produces  odd  symmetry  in  the  field 
about  the  pressure  release  boundary  which  must  now  satisfy 

'P(-z  +  2;7(x,y),x,y)  =  -'P(z,x,y).  (5) 
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Thus,  the  image  oeean  must  satisfy  the  SPE  of  the  form 
ST*  d  T)  ST* 

^  ^  ^  i'^iDSPE  ~^^3DSPE^^  ’  (6) 

OX  ox  oz 

Tappert  andNghiem-Phu  (1985)  introdueed  a  field  transformation 


'P(z,x,y)  =  <^ 


^(z,x,y) 


^{z,x,y) 


,  z  >  r]{x,y) 
z  <  ?7(x,y). 


(V) 


In  the  real  ocean,  this  transformation  produces 


a'P 


^  =  -K(^3DSPE  +  U.oSPeW  ’  2 


(8) 


In  the  image  ocean,  we  must  first  evaluate  the  left-hand  side  of  that  image  ocean 
form  mentioned  previously  in  Eq.  (6) 


S'P 


dx 


dx  dz 


=  e 


dx 


d^Tj 


- ilk^ - :r{z-Tj)^'  -i2k^ 

ox  ox 


m  -5;7  5'P 


ydx  j 


^  +  2- 


dx  dz 


(9) 


so 


The  right-hand  side  of  Eq.  (6)  includes  the  operator 

1  a'  1  5" 


3DSPE 


2k^  dz^  2e  dy 

o  o 


2  a,.2  ’ 


(10) 


1  — (z-17) 

T  W  = _ p  Sx 

3DSPE  ^  ^  ,  2 

2k 


dx  dz 


drj 

ydx  j 


'P 


+  - 


e^'p 


-i4k„ 


^  d^rj  ^  ^  drj  drj 

- -{z-rf) - 

y  dxdy  dx  5y 


dy 


-i2k„ 


(r  7])  2  i2k 


dxdy 


dxdy  dy  dx  d  y 


^  d^n  ,  ^  drj  drj 

- -{z-Tj) - 

dxd)y  dx  d2y  j 


2  A 


'P 


(11) 
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After  combining  the  right-hand  side  with  the  left-hand  side,  we  get 


dx 


=  -ik,. 


1  1 

( a^'F 

a"'F^ 

^d^n  . 

1 

V 

+ 

^3DSPE 

d^rj 

dxdy 


{Z-T])- 


dij  dij 
dx  dy 


a'F 

dy 


dxdy^  dxdy  dy 


drj  d^rj 
dx  dy^ 


-i2k„ 


''  d^rj 
^dxdy 


{z-rt) 


drj  dri^ 
dx  dy  j 


'F  . 


(12) 


Consequently,  we  can  combine  these  expressions  together  into  a  single 
formulation, 

a'F  ~  ~  ~ 

- ^-ik{T  +f/  )'F  (13) 

Qy^°  °P°P  ^  ' 

with 


and 


f  = 

op 


2DSPE 


2DSPE 


a 

+ — 
k 


d^ri 

dxdy 


{Z-T]) 


drj  drj  d 
dx  dy  dy 


,z>  r/(x,y) 
,  z  <  r/(x,y) 


(14) 


^op=< 


U. 


3DSPE 


,  z>  ri{x,y) 


U,dspe-^^{^-J1) 


(15) 


-F- 


d^t]  d^t]  dt]  dt]  d^t] 

-{z-rj)-! - 


k  d^rj 


dxdy^ 


dxdy  dy  dx  dy 


dij  dr] 


(z-7)-- 
dxdy  dx  dy 


z  <  ;7(x,y) 
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The  previously  mentioned  equations  represent  an  exaet  formulation  of  the  3-D 
rough  surfaee  seattering  utilizing  the  SPE  approximation,  as  originally  derived  by  Smith 
(2012).  The  form  of  the  3-D  kinetie  energy  operator  introduces  some  computational 
complexities.  Specifically,  the  field  above  the  interface  undergoes  a  separate  cross-range 
operator.  And  while  it  is  possible  to  apply  this  operator,  it  is  not  possible  to  apply  it  at  a 
single  depth  mesh  point  since  the  rough  surface  varies  with  cross-range  (i.e.,  ri(x,ym)  7^ 
Tl(x,yn)  )• 

In  order  to  overcome  this  difficulty,  an  approach  has  been  suggested  by  Smith 
(2013).  First,  at  each  range  step  (xi),  the  3-D  field  is  transformed  from  cross-range  spatial 
variable  (y)  to  cross-range  wavenumber  variable  (ky)  over  all  depths.  Concurrently, 
transformation  from  depth  variable  (z)  to  vertical  wavenumber  variable  (kz)  should  also 
be  performed.  Second,  a  copy  of  the  field  in  (ky,kz)  space  should  be  implemented.  Third, 
to  one  copy  of  the  field,  the  for  z  >  p  should  be  applied,  while  to  the  other  copy  the 

for  z  <  p  should  be  applied.  Fourth,  each  copy  is  transformed  back  to  the  spatial  (y,  z) 

domain.  Fifth,  based  on  the  local  surface  height  at  range  (xi),  and  for  each  cross-range 
(ym),  the  appropriate  form  of  the  field  should  be  recombined.  Finally,  the  appropriate 
form  of  should  then  be  applied  to  complete  the  range  step  to  (xi+i). 

The  preceding  3-D  expressions  were  developed  using  the  SPE  approximation  in 
Cartesian  coordinates  (x,  y,  z).  The  development  of  the  equations  in  cylindrical 
coordinates  (r,  0,  (j))  closely  follows  the  previous  approach,  beginning  with 

=  ,z>n(r4),  (16) 

where 


T  =T 

^op  ^1,DSPE 


1 


1 


Ik^  dz^  2ky  dcj) 


,U  =U. 

^  op  . 


3DSPE  2 


■1). 


(17) 
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This  time  the  rough  surface  boundary  condition  is  introduced  by  defining  the 
location  of  the  pressure  release  rough  surface  at  depth 

z  =  ri{r,(f)).  (18) 

Then  the  pressure  release  boundary  condition  is  defined  by 

'P(z  =  ;7(r,^/i),r,(^)  =  0.  (19) 


Introducing  an  image  field  for  z  <  ;/  produces  the  symmetry  defined  by 

'P  (-Z  +  Irjir,  (j)\  r,  (f))  =  -'P  (z,  r,  (f)).  (20) 

Thus,  the  image  ocean  must  satisfy  a  PE  of  the  form 

f^P  d  t]  fT* 

+  2^—  =  -ikST^^+UJ^>  ,z<;7(r,^/i).  (21) 

or  or  oz 

As  before,  we  introduce  a  field  transformation  of  the  form 


'P(z,r,^zi) 


'P(z,r,^) 


,  z>  ijir,^) 
,  z  <  ;7(r,^/i). 


(22) 


12 


These  equations  follow  exaetly  the  same  development  as  before,  but  now  we  end 
up  with  the  eombined  formulation 


with 


a'T 

- ^-ik{T  +U  )'T 

n  ^  on  on ' 

or 


o  ^  op 


op' 


and 


To,=< 


hDSPE 


T  + 

^3DSPE  ^ 


a 


k„r 


d^n  ^  ^  dn  dn 

- -(z-n) - 

drd^  dr  d(l) 


A 

dcj) 


,z>  77(r,^^) 
,  z  <  77(r,^/i) 


(23) 


(24) 


^3DSPE  ’  ^  ^ 


u. 


3DSPE 


+  - 


kr^ 


d  i]  d  ij  dij  dij  d  ij 

-  (z  -  77)  -  2 - 7  -  i2k_ 


(25) 


drd(^ 


drd^  5^  dr 


dr]  drj  drj 

- (z-7) - 

\^drd(l)  dr  drj) j 


,z<  ri{r,^) 


The  implementation  of  these  equations  follows  the  same  outline  as  provided  in  the 
Cartesian  coordinates  approach. 

However,  it  has  been  shown  (Smith  and  Colosi  2008)  that  the  SPE  approximation 
performs  poorly  in  3-D  coordinates  since  cross-ranges  separated  from  the  center 
introduce  lateral  angles  that  are  not  well  treated.  They  suggest  that  a  better  approach  in  3- 
D  is  to  utilize  the  2-D  wide-angle  parabolic  equation  (WAPE)  approximation  in  order  to 
produce  the  proper  3-D  structure  of  the  field.  It  is  shown  that  several  problems  in  sound 
propagation  cannot  be  studied  by  using  a  small-angle  parabolic  equation,  and  need  a 
WAPE  for  their  solution  (Ostashev  and  Blanc-Benon  1997). 
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Smith  (2012)  discussed  the  approach  to  these  scattering  corrections  using  a  wide- 
angle  PE.  He  noted  that  the  same  seattering  eorreetions  result  if  third  order  derivatives  of 
the  field  and  surfaee  slope  are  ignored.  Therefore,  we  shall  utilize  the  3-D  WAPE 
propagator  funetions  with  the  seattering  eorreetion  terms  derived  previously.  These 
operators  are  defined  in  cylindrical  coordinates  by 


Meanwhile,  the  WAPE  introduees  new  terms  whieh  depend  upon  vertical 
derivatives  of  the  field  of  seeond  order  and  higher  and/or  produets  with  terms  that  are 
second  order  and  higher  factors  of  the  surface  slope.  Even  though  the  WAPE 
implementation  in  the  water  column  allows  for  better  accuracy  of  high  angle  propagation, 
the  existing  application  of  the  rough  surface  scattering  approach  does  not  appear  to 
support  the  same  level  of  high  angle  seattering  and  assumes  small  surfaee  slopes. 

An  approximation  is  done  to  the  in  the  image  ocean  by  dropping  the  third 

order  and  higher  terms  in  the  derivative  of  the  surface  displacement  that  negleets  the 
effeets  of  very  high  slopes  in  the  scattering  and  the  third  order  and  higher  azimuthal 
derivatives  that  invoke  a  small-angle  approximation.  However,  seeond  order  derivative 
terms  are  maintained.  This  is  eonsistent  with  seeond  order  derivatives  of  the  propagating 
field,  although  sueh  terms  are  not  diffieult  to  include. 
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III.  NUMERICAL  IMPLEMENTATION 


A,  HYBRID  IMPLEMENTATION  USING  A  TRI-DIAGONAL  SOLUTION 
FOR  THE  SCATTERING  TERM 

The  forms  of  the  derived  equations  using  cylindrical  or  Cartesian  coordinates  are 
identical  with  the  obvious  substitution 


d  Id 

- ^ - . 

dy  r  d(p 


(27) 


For  conciseness,  we  shall  only  consider  the  Cartesian  form  in  what  follows, 
although  the  model  implementation  is  done  using  cylindrical  coordinates. 

The  derived  equations  in  the  previous  chapter  contained  the  additional  term  in  the 
r^operator  for  z<  7(^x,yjwhich  cannot  be  applied  in  the  Ar^wavenumber  domain,  as  is 

traditionally  done  for  .  Specifically,  if  we  write 


T  ={ 

Op  * 


'-?>DWAPE 


dy 


,z>  r]{x,y) 
,  z<  ri{x,y)  , 


(28) 


we  see  that  =  T^dwape  ( — 7^ — 7)’  which  can  be  applied  in  the  (ky,kz)  wavenumber 

dy  dz 

d  d 

domain.  However,  7? —  =  R(x,y,z) —  cannot  be  treated  so  easily. 
dy  dy 

d  d  d 

It  is  the  combined  dependence  of  7? —  on  both  y  and  — (or  ^and  — in  the 

dy  dy  d(p 

cylindrical  coordinates)  that  does  not  allow  a  simple  implementation  via  a  Split-Step 
Fourier  (SSF)  algorithm.  Instead,  the  following  hybrid  approach  is  proposed. 
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Loosely  following  the  approach  of  Yevick  and  Thomson  (1997),  we  will  focus  on 
the  image  field  evolution  equation 


dW  ~  ~ 
—  ^-ikH  'T 
ax  ° 


where,  is  a  Hamiltonian-like  operator  that  defines  the  evolution  of  the  PE  field 
function  in  range.  We  may  then  use  the  approximate  marching  algorithm  solution 


(29) 


^^(x-l-  Ax)  =  e 


-ik  ArH  (r) 


'P(x). 


(30) 


Here, 

.  (31) 

so  we  may  define, 

„  -ik^AxR—  ^ 

^(x  +  Ax)  =  e  ^2) 

Note  that  this  form  assumes  that  the  rough  surface  correction  is  applied  after  the 
kinetic  energy  operator  (r3^). 


The  Monterey-Miami  Parabolic  Equation  (MMPE)  model  employs  a  centered- 
step  scheme  implementation,  which  in  the  real  ocean  takes  the  form  (Smith  and  Tappert 
1993,  Smith  2001) 


'P(x  +  Ax)  =  e  ^  FFTy^  <e 


'FFT... 


-i—k^AxU-tn 

e  2°  ^>(x) 


(33) 


where  'P(x)  = 'P(x,y,z)and  the  scalar  form  of  the  operator  the  (ky, 

kz)  wavenumber  space. 
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The  image  ocean  solution  can  utilize  the  SSF  algorithm  for  the  term,  but  must 


employ  an  alternative  method  for  the  R — term.  Thus,  we  will  implement  a  centered-step 

dy 

scheme  defined  by 


'T(x  + Ax)  =  e 


—  1  I  “  Ax  7^  D  ( ^  V  ) 
y2  -  -  -yz 


2  dy ppjz-i 


-i—k^^AxU^j^  (x) 


T(x) 


>.  (34) 


In  order  to  implement  this  equation,  we  should  determine  how  to  apply  the 


operator  e  ^  with  g  =  k  Ax  and  y  -  —  by  defining  the  intermediate  step  solution  as 

"  dy 


'T'(x  +  Ax')  =  FFT-'  I e-‘^^^°FFT^^ 


?  2  ’“'V(x) 


(35) 


The  final  solution  at  the  end  of  each  range  step  is  then  obtained  from 


-N-  -i—cU-iry(x+Ax)  ^ 

'P(x  +  Ax)  =  e  2  'PCx  +  Ax'), 


(36) 


where 


W(x  +  Ax')  =  e-‘^‘'’'W'(x  +  Ax'). 


(37) 
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The  rough  surface  correction  operator  can  be  represented  by  the  Fade 
approximant  (Collins  1990) 


(38) 


which  is  valid  for  small  y ,  implying  small  corrections  in  cross-range  scattering  at  each 
step.  From  this  approximation,  we  obtain  a  finite  difference  system  of  equations  based 
upon 


1 

1 

II 

-t 

'F'(^  +  Ajr')- 


This  may  be  written  in  a  vector-matrix  form  as 


A^  =  B^'  =  C  , 


where 


^2 

'i'3 

'F'  = 

'f; 

'i'., 

'i".. 

N  ' 

N 

and  W^=W(yJ  ,  n=l,2,3....N. 


(39) 


(40) 


(41) 
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The  expressions  required  for  the  ^4  and  ^matriees’  components  will  be  based  on  a 
3 -point  centered  approximation  to  the  derivative 

~  1  ~  ~ 

/!'  =  —  (42) 

dy  2Ay 

^  0 

where  =  'F(jv  =  T  -  Ah)  -  the  end  points  (n=l  or  n=N),  we  shall  invoke  the  2-point 


approximations 


~  S'F  1  ~  ~ 

/i'  =  — 

dy  Ay  “ 


~  S'F  1  ~  ~ 

/F  = -  ('F-'FJ. 

0y  Ay'  ' 


The  former  is  used  if  n=N  and  the  latter  applies  when  n=l . 
The  matrices  may  be  defined  in  terms  of  a  and  b  where 


\  +  UcRy  T  ='F  +-icRy^  = 

2^'  o  O  ^2^  ^  •  o 


'F  (T  ,)  =  'F  +^('F 

°  2  2Ay  *  °  4Ay  ' 


from  which  we  can  introduce 


1 

a  = - - 

"  4  Ay 


,  -1 

b  = - ^ 

"  4  Ay 


R  =  R{y  ). 

n  n-' 
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Consequently,  we  ean  define  two  different  matrices 


and 


"  (l-2aj)  2aj  0  0  0  0  ^ 

-a^  1  0  0  0 

0  -^3  1  ^3  0  0 

0  0  0  0  0 

0  0  0  -a^_^  1 

^  0  0  0  0  -2a^  (l  +  2a^)  ^ 


"  (l-2Z>j)  2Z>j  0  0  0  0  ^ 

1  Z>2  0  0  0 

0  -Z>3  1  Z>3  0  0 

0  0  0  0  0 

0  0  0  -K-.  1  K-r 

,  0  0  0  0  -2b^  (1  +  26J^ 


Now  the  matrix  equation  becomes 


A^>  =  C  = 


c, 

C3 


(47) 


(48) 


(49) 


C 


N 


with  the  elements  of  C  are  computed  from  .  Equation  (50)  can  be  solved  by  a  tri¬ 
diagonal  matrix  algorithm  (Press  et  al.  2007)  according  to  the  recursive  relations 


N 


c 

'-N 


'P 


C'-a^ 


M 


1=  1,2,3, 


,N-2,N-1  , 


(50) 
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where 


2a^ 

l-2aj 

a. 

_ l _ 

1  +  a'  ,a. 

l-l  l 


i=l 


i=2,3,...,N-l 


l-2aj 


C.  +  C  .Q. 

1 _ i-\  I 

\  +  a'  ,a. 

l-\  l 


i=l 


i=2,3,...,N-l 


C  -\-2ci  C 

_ AT _ N  N-\ 

(l  +  2a^)  +  2a^a;_j 


i=N 


(51) 


(52) 


The  outeome  of  this  tri-diagonal  matrix  solution  is  then^(x  +  Ax')  to  whieh  we 

(x+Ax) 

simply  apply  the  sealar  veetor  multiplier  e  ^  to  obtain  the  solution  at  the  end  of 

the  range  step  'T(x  -I-  Ax) . 
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B, 


DEFINITION  OF  ROUGH  SURFACE  PERTURBATIONS 


In  order  to  appropriately  test  the  3-D  scattering  model  just  developed,  it  is 
required  that  we  analytically  define  some  simple  surface  perturbations.  In  this  section,  we 
shall  define  five  simple  surface  perturbations.  These  will  be  defined  in  cylindrical 
coordinates,  since  the  model  is  implemented  this  way. 

First,  we  define  a  surface  with  constant  slope  in  order  to  determine  the  way  the 
different  radials  will  be  affected  by  that  surface  slope  using  the  equation 


rj {r,  (f))  =  20 +  0.l{r  cos  (ft),  (53) 

where  T]{r,(j))  is  the  surface  height  with  respect  to  range,  r  ,  and  the  expressions  for  the 
required  derivatives  are  then 


—  =  O.lcos(^zJ), 
dr 

^  =  0 

dr^ 

|^  =  -0.1rsin((zi). 


d^rj 

drd^ 


=  -O.lsin(^zi). 


(54) 


This  surface  is  depicted  in  Fig.  1. 

Second,  we  define  the  surface  roughness  for  concentric  sinusoidal  waves 

?7(r,  (p)  =  G  cos{kr),  (55) 


where  //(r,^)  is  now  the  surface  wave  height  with  respect  to  range,  r  (independent  of 
bearing  (j)  in  this  case),  cr  is  the  surface  wave  height  amplitude,  and  k  is  the  surface 
wavenumber  {luj  where  is  the  wavelength  of  the  surface  perturbation).  Notice  here 

that  since  the  surface  does  not  depend  upon  the  bearing  angle,  the  propagation  out  along 
all  radials  should  be  affected  similarly.  Then 
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drj 

dr 

d^Tj 


=  -k(j  sm{kr). 


=  -k^a  cos{kr), 


^  =  0, 

d(p 


drd^ 


This  surface  is  depicted  in  Fig.  2. 


(56) 
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Figure  2.  Concentric  sinusoidal  surface  perturbations. 


Third, 

unidirectional 

where  (j)  is  the 


we  define  the  equations  introducing  the  surface  roughness  for 
sinusoidal  waves  by 


T]{r,  ^')  =  <7  cos{kr  cos  (p), 

angle  of  propagation  along  each  radial.  Then 

drj 
dr 
d^rj 


(57) 


=  -A:cr  cos(^zi)  sin[A:rcos(^zi)], 

=  -(7{k  cos(^zi)]^  cos[A:r  cos(^zi)], 


drj 

—  =  kro  sin((zi)  sin[A:r  cos((zi)], 
d(j) 


(58) 


d^ri 

drd(l) 


=  k(j  sin(^zi)  {sin[^  cos(^zi)]  -l-  kr  cos(^)  cos[A:r  cos(^zi)]} . 
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Figure  3.  Planar  sinusoidal  surface  perturbations. 


Fourth,  we  define  a  single  upward  sinusoidal  perturbation  on  one  side  only  by 
using  the  same  surface  roughness  equation  as  the  unidirectional  surface  waves,  but  this 
time  we  limit  the  angle  of  propagation  (j)  to  be  within  the  forward  region  of  the  source  and 
only  allow  one  oscillation, 

-90  <^<9Q, 

r](r,^)  =  <jcos{krcosg))  ,  (59) 

r  cos 
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Elsewhere,  the  surfaee  will  be  eonstant  with  a  value  for  rjir,^')  equal  to  the 
maximum  surfaee  wave  height  cr.  The  derivatives  in  the  region  of  the  displaeement  are 
then 


—  =  -ka  eos({z5)  sm[  A:r  eos(^zi)] , 
dr 


=  -(j{k  cos(^zi)]^  cos[A:r  cos(^zJ)], 


d^Tj 

—  =  A:rcrsm(^zi)sm[A:r  eos({z))], 

d^ 


(60) 


d^rj 

drd^ 


=  k(j  sm(^z))  {sm[  A:r  cos(^zJ)]  +  kr  cos{(/))  eos[  A:r  eos(^z5)]} , 


while  over  the  other  direetions  and  ranges,  the  sea  surfaee  is  defined  as 

Tj{r,(j))  =  a, 

^  =  0, 

dr 


d^rj 

dr^ 


=  0, 


(61) 


^=0, 

d^ 

d^rj 


drd<f) 


=  0. 


This  surfaee  is  depleted  in  Fig.  4. 


Finally,  we  flip  the  single  upward  sinusoidal  perturbation  wave  and  shift  it  until  it 
reaehes  the  mean  sea  level,  where  //(r,  ^)  equals  zero  everywhere  exeept  in  the  forward 
region  of  the  souree,  where  the  sea  surfaee  equation  is  defined  by 


-90  <  <  90, 

ri{r,<f))  =  G\V-QO%{krQO%<f))\  ,  (62) 

rcos(p<A^. 
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Figure  4.  Upward  sinusoidal  surface  perturbations. 


The  derivatives  in  the  region  of  the  displacement  are  then 


dvi 

— -  =  kcr  cos{(f>)sm[kr  cos{(f>)\, 
dr 

=  (j\k  cos(^/i)f  cos[^r  cos(^)], 
dr 


— -  =  -krcT  sin(^)  sin[^r  cos(^)] , 
d(j) 


^  =  -k<j  sin(^/i)  |sin[^r  cos(^/i)]  +  kr  cos(^)  cos[kr  cos(^)]} . 


while  all  over  the  other  directions,  the  sea  surface  is  defined  as 


(63) 
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Tj{r,(j))  =  ^, 


dij 

dr 


=  0, 


d^ri 

dr^ 


=  0, 


drj 

d(j) 


=  0, 


drd(j) 

This  surface  is  depicted  in  Fig.  5. 


(64) 


Figure  5.  Downward  sinusoidal  surface  perturbations. 
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IV.  MODEL  RESULTS 


A,  CONSTANT  SLOPE  SURFACE  VERSUS  CONSTANT  SLOPE  BOTTOM 

Although  a  sloping  surface  is  a  relatively  benign  definition  of  a  perturbed  surfaee, 
it  does  allow  a  eomparison  of  the  rough  surfaee  algorithm  defined  here  with  a  formally 
equivalent  environment  defined  by  a  sloping  bottom.  The  solution  to  a  eonstant  slope 
surfaee  is  identieal  to  a  eonstant  slope  bottom  through  a  simple  environmental  rotation, 
assuming  ideal  point  sourees  in  both  oases.  The  approach  to  the  3-D  problem  with  a  non- 
oonstant  bottom  is  no  different  from  the  approaoh  to  the  2-D  problem,  where  the  MMPE 
model  simply  adjusts  the  bohom  depth  at  eaoh  range  step  and  radial. 

1.  Constant  Slope  Surface  and  Flat  Bottom 

Solutions  to  the  problem  of  a  sloping  surfaee  over  a  fiat  bottom  are  depleted  in  the 
following  figures.  The  souroe  is  looated  at  a  depth  of  50  m  and  the  bottom  depth  is 
highlighted  with  a  white  line  at  120  m  depth.  Figures  6  and  7  display  the  solutions  along 
opposing  radials,  with  the  former  along  the  direetion  of  deseending  surfaee  depth  and  the 
latter  along  the  direetion  of  aseending  surfaee  depth.  In  both  oases,  the  bottom  interfaoe  is 
observed  to  be  flat. 

In  Figure  8,  the  transmission  loss  at  a  single  depth,  oorresponding  to  the  souroe 
depth,  is  extraoted  from  all  radials.  Note  that  in  the  ease  of  a  perfeotly  fiat  surfaee  and 
bottom,  this  would  produoe  oonoentrio  rings  of  interferenoe  patterns  in  the  field.  In  this 
ease,  however,  the  direotionality  of  the  propagation  along  different  radials  is  apparent. 
Transforming  this  data  to  a  polar  ooordinate  system  produoes  the  results  presented  in 
Figure  9.  The  effeot  observed  is  oonsistent  with  general  expeotations  of  a  sloping  surfaee, 
where  the  propagation  in  the  direetion  of  the  deseending  surfaee  produoes  interferenoe 
patterns  oompressed  in  range  and,  alternatively,  the  propagation  in  the  direetion  of  the 
aseending  surfaee  produoes  interferenoe  patterns  elongated  in  range. 
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Figure  6.  0-degree  radial  for  a  source  at  depth  of  50  m. 


Transmission  Loss  (dB  re  1m) 


Range  (km) 

Figure  7.  180-degree  radial  for  a  source  at  depth  of  50  m. 
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Transmission  Loss  (dB  re  1m) 


Radial  (deg) 

Figure  8.  Transmission  loss  pattern  at  source  depth. 


Range  (Km) 

Figure  9.  Constant  slope  surface  effect  on  the  waves  emitted  by  a  50  m  deep  source. 
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2,  Constant  Slope  Bottom  and  a  Flat  Surface 

The  corresponding  solution  is  now  computed  by  interchanging  the  surface  and 
bottom  characteristics.  The  environment  now  has  a  smooth  flat  surface  with  a  sloped 
bottom  which  should  produce  the  equivalent  field  solution,  albeit  rotated,  as  in  the 
constant  slope  surface  case.  In  this  case,  the  surface  depth  is  set  to  a  constant  value  of 
zero  and  the  source  depth  is  shifted  upwards  to  30  m  depth  (consistent  with  source  depth 
relative  to  the  surface  in  the  previous  case),  while  the  average  depth  of  the  medium  will 
be  100  m.  The  2-D  bottom  bathymetry  for  this  environment  is  defined  as 

bdint  =  bdint-0.1(rcos^zi)  (65) 

where  “bdint”  is  the  interpolated  bottom  depth  at  each  range  step  and  (j)  is  the  bearing 
angle  of  propagation. 

Figures  10  and  11  present  the  results  along  opposing  radials.  The  results  in  these 
figures  should  compare  directly,  through  a  simple  rotational  transformation,  to  the  results 
presented  earlier  in  Figures  6  and  7. 


Transmission  Loss  (dB  re  1m) 
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Figure  10.  0-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  11.  1 80-degree  radial  for  a  source  at  depth  of  30  m. 

3,  Constant  Slope  Surface  Model  versus  Constant  Slope  Bottom  Model 

In  order  to  compare  results  from  the  two  different  approaches,  a  common  depth 
for  extracting  transmission  loss  must  be  defined.  The  most  convenient  depth  to  choose  in 
both  cases  is  the  depth  of  the  bottom  bathymetry.  Other  than  some  geometrical  rescaling 
(due  to  rotation)  of  the  range  from  the  source,  the  transmission  loss  traces  along  the 
bottom  depth  should  agree  between  these  two  methods. 

Figure  12  shows  the  initial  results  directly  extracted  from  each  solution  at  the 
bottom  interface  in  polar  coordinates.  Due  to  the  inclination  presented  by  the  sloping 
bottom  model,  we  must  adjust  the  range  defined  in  the  calculation  in  order  to  compare 
results.  Figures  13  and  14  utilize  this  range-scaling  correction  to  compare  the 
transmission  loss  traces  along  opposing  radials. 

As  can  be  noted  in  Figures  13  and  14,  the  comparison  between  solutions 
computed  with  a  2-D  sloping  bottom  using  the  standard  3-D  MMPE  and  this  new 
algorithm  with  a  2-D  sloping  surface  show  good  agreement. 
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Figure  12. 


Model  solutions  for  the  water-bottom  interfaee:  (left)  the  constant  slope  surface; 

(right)  the  constant  slope  bottom. 


Figure  13.  Transmission  loss  at  water-bottom  interface  along  the  0  radial. 
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30 


Range(Km) 

Figure  14.  Transmission  loss  at  water-bottom  interfaee  along  the  180  radial. 

B.  CONCENTRIC  SINUSOIDAL  SURFACE  PERTURBATIONS 

As  an  alternative  approaeh  to  eonfirm  proper  model  behavior,  we  now  examine 
the  results  derived  from  a  surfaee  perturbation  defined  as  concentric  sinusoidal  waves. 
Theoretically,  all  the  radials  should  be  affected  in  the  same  manner  due  to  the  defined, 
radially-independent  environment. 

Figures  15  and  16  present  the  transmission  loss  field  results  from  opposing 
radials.  As  expected,  the  results  appear  identical  in  their  structure  along  each  radial. 
Figure  17  displays  the  transmission  loss  results  extracted  from  a  constant  depth  of  30  m, 
the  same  as  the  source  depth.  The  results  transformed  to  polar  coordinates  are  presented 
in  Figure  18. 

The  results  show,  as  expected,  that  the  acoustic  field  will  propagate  identically 
along  all  bearing  angles,  producing  the  circular  patterns  shown  in  Figure  18. 
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Figure  15.  0-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  16.  1 80-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  17.  Transmission  loss  pattern  at  souree  depth. 
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Figure  18.  Concentric  sinusoidal  surface  perturbations  effect  on  a  transmission  loss  field 

extracted  at  source  depth. 
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c. 


PLANAR  SINUSOIDAL  SURFACE  PERTURBATIONS 


We  now  examine  the  results  computing  using  the  planar  sinusoidal  surface 
perturbations.  While  the  forward  and  backward  (0  and  180  degrees)  radials  should 
produce  the  same  results  due  to  symmetry,  there  should  be  significant  differences 
between  those  and  other  radials,  where  the  surface  perturbation  along  the  path  is 
different. 
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Figure  19.  0-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  20.  1 80-degree  radial  for  a  source  at  depth  of  30  m. 
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Figures  19  and  20  display  the  symmetric  radials  along  0  and  180  degrees.  As 
expected,  the  transmission  loss  patterns  appear  identical  between  the  two  bearings. 
However,  since  the  propagation  will  respond  differently  along  each  radial,  an 
examination  of  the  field  at  a  single  depth  produces  the  results  displayed  in  Figures  21 
and  22. 

Note  that  the  dominant  scattering  features  present  (most  notably  in  Figure  22) 
occur  upon  the  first  interaction  with  the  downward  perturbation  of  the  surface  roughness. 
The  impact  on  the  field  is  clearly  parallel  to  the  plane  of  the  surface  perturbation,  and  for 
ranges  shorter  than  +/-  0.1  km  on  the  horizontal  cross-range  axis,  energy  appears  to  be 
channeled  between  the  downward  portions  of  the  surface  roughness  on  either  side. 

Beyond  +/-  0.1  km  on  the  horizontal  cross-range  axis,  the  impact  of  the  surface 
roughness  is  less  noticeable.  This  is  due  to  the  dominant  scattering  occurring  within  the 
first  +/-  0.1  km.  Beyond  this  range,  much  of  the  surface  interaction  has  been  shadowed  by 
the  near  field  scattering. 
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Figure  2 1 .  Transmission  loss  pattern  at  source  depth. 
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Figure  22.  Planar  sinusoidal  surface  perturbations  effect  on  a  transmission  loss  field 

extracted  at  source  depth. 

D,  UPWARD  SINUSOIDAL  SURFACE  PERTURBATIONS 

In  order  to  isolate  some  of  the  scattering  phenomena  occurring  from  the  previous 
type  of  planar  sinusoidal  surface  perturbation,  and  to  break  the  symmetry  of  the  problem 
to  ensure  proper  model  behavior,  the  next  surface  to  test  is  simply  a  single  sinusoidal 
surface  displacement.  We  start  initially  by  creating  a  single  upward  sinusoidal  surface 
perturbation  to  the  field,  identical  in  form  to  the  previous  surface  on  one  side  only  and 
out  to  0.1  km.  We  select  the  0  radial  to  represent  the  source  front-side  of  the  acoustic 
propagation  and  the  180  radial  to  show  the  response  on  the  back-side  of  the  source. 

The  results  are  presented  in  Figures  23  through  26.  As  expected,  the  back-side  of 
the  propagation  shows  little  impact  from  the  surface  perturbation  on  the  front-side.  There 
may  be  some  minor  effect  along  the  boundary  separating  front-side  and  back-side,  but 
this  is  to  be  expected  due  to  horizontal  scattering.  In  addition,  the  scattering  in  the  front¬ 
side  direction  is  very  similar  to  the  previous  multiple  sinusoidal  planar  perturbation. 
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confirming  that  most  of  the  scattering  effects  are  due  to  the  initial  interactions  with  the 
downward  sloping  portion  of  the  surface. 


Figure  23 .  0-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  24.  1 80-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  25.  Transmission  loss  pattern  at  source  depth  of  30  m. 
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Figure  26.  Upward  sinusoidal  surface  perturbations  effect  on  a  transmission  loss  field 

extracted  at  source  depth. 
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E,  DOWNWARD  SINUSOIDAL  SURFACE  PERTURBATIONS 

As  an  alternative  to  the  single  sinusoidal  planar  perturbation  presented  previously, 
we  also  ereated  a  similar  surfaee  that  perturbs  the  surfaee  downward  initially  before 
returning  upward  to  a  eonstant  depth.  One  of  the  primary  differences  of  this  perturbation 
was  that  it  places  the  strongly  scattering  downward  slope  of  the  surface  perturbation 
immediately  adjacent  to  the  source  range.  This  would,  therefore,  be  expected  to  introduce 
stronger  scattering  at  higher  angles,  particularly  in  the  near-field.  As  before,  we  choose 
the  0  radial  to  represent  the  front-side  path  of  the  acoustic  propagation  and  the  180  radial 
to  represent  the  back-side  field  of  the  source. 

It  is  worth  noting  that  during  early  testing  of  the  model’s  stability,  this  particular 
surface  perturbation  initially  generated  non-physical  results.  Closer  inspection  revealed 
that  the  earlier  version  of  the  model  was  utilizing  an  SPE  approximation  to  compute  the 
azimuthal  wavenumber  propagator  function,  resulting  in  large  scattering  angles  being 
misrepresented  in  the  solution.  As  discussed  in  the  theory  section  of  this  thesis,  it  is  quite 
important  to  employ  wide-angle  approximations,  even  in  cross-range/bearing-angle,  in 
order  to  properly  compute  the  azimuthal  field  structure.  After  upgrading  the  code  to 
account  for  this,  the  results  from  this  perturbation  model  became  physically  realizable. 
(All  results  reported  in  this  thesis,  therefore,  utilized  the  WAPE  approximation  for  all 
operators.) 

Eigures  27  through  30  display  the  computed  transmission  loss  field  results  for  this 
surface  perturbation.  Again,  the  back-side  of  the  field  shows  little  impact  from  the 
perturbation,  except  very  near  the  edge  separating  the  front-side  from  the  back-side, 
which  is  to  be  expected.  In  addition,  the  front-side  scattering  is  much  stronger  in  this  case 
than  in  the  previous  one.  There  is  also  clearly  diffraction  of  the  field  filling  in  the 
classical  shadow  zone  on  the  back-side  of  the  downward  surface  perturbation. 

It  is  natural  to  expect  that  the  scattering  from  the  front-side  would  eventually  be 
transferred  across  the  boundary  to  the  back-side  of  the  field.  Because  of  the  limitation  on 
the  source  angle  in  the  vertical,  and  the  relatively  short  range  of  the  calculations,  this  was 


43 


Depth  (m)  Depth  (m) 


not  clearly  observed  in  any  of  these  calculations.  It  is  worth  exploring  this  issue  further  in 
future  analysis. 


Transmission  Loss  (dB  re  1m) 


Range  (km) 

Figure  27.  0-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  28.  1 80-degree  radial  for  a  source  at  depth  of  30  m. 
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Figure  29.  Transmission  loss  pattern  at  souree  depth. 
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Figure  30.  Downward  sinusoidal  surfaee  perturbations  effect  on  a  transmission  loss  field 

extracted  at  source  depth. 
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F. 


COMPARISON  OF  THE  EFFECT  OF  3-D  VERSUS  N  x  2-D 


One  of  the  goals  of  this  projeet  is  to  develop  an  aeeurate  means  of  computing  the 
3-D  scattering  from  a  rough  2-D  surface,  and  then  comparing  solutions  between  this  3-D 
model  and  standard  Nx2-D  models.  Such  Nx2-D  models  are  significantly  more  efficient 
and  require  much  less  in  the  way  of  computing  power.  And  so,  understanding  when  a  true 
3-D  solution  is  required  will  help  avoid  unnecessary  computational  complexity.  In  this 
section,  we  will  compare  results  computed  between  an  Nx2-D  version  of  the  MMPE 
model  that  utilizes  a  consistent  rough  surface  scattering  algorithm  with  the  3-D  version 
described  in  this  thesis. 

1,  Constant  Slope  Surface  Perturbations 

Utilizing  the  same  environment  defined  previously  for  the  constant  sloping 
surface,  the  results  from  the  3-D  model  and  the  corresponding  Nx2-D  version  were 
computed  out  to  a  maximum  range  of  0.6  km.  The  goal  of  carrying  out  these  calculations 
to  a  slightly  longer  range  than  before  is  simply  to  allow  for  more  azimuthal  coupling  to 
occur  which  would  differentiate  the  impact  of  the  3-D  calculation.  Figures  31  through  34 
present  the  model  results  from  this  environment  along  the  different  radials  and  at  a 
constant  depth. 

The  results  from  the  two  models,  particularly  as  displayed  in  Figure  34  in  polar 
coordinates,  show  little  difference  between  the  models.  This  is  perhaps  not  too  surprising, 
since  azimuthal  coupling  is  known  to  have  a  larger  impact  at  longer  ranges.  More 
significant  differences  would  likely  be  observed  if  the  calculations  were  carried  out  at  a 
lower  frequency  and  the  maximum  range  was  increased  significantly.  However,  due  to 
the  nature  of  the  sloping  surface,  longer  ranges  would  cause  the  surface  to  cross  the 
bottom  interface,  generating  nonsensical  environments. 

In  an  effort  to  highlight  the  differences  computed  from  the  3-D  and  Nx2-D 
models.  Figure  35  displays  the  difference  in  transmission  loss  values  at  a  constant  depth. 
As  noted  previously,  the  differences  are  relatively  minor,  and  only  begin  to  become 
prominent  at  longer  ranges.  Furthermore,  the  strongest  differences  occur  along  the  radials 
at  90  and  270  degrees,  corresponding  to  bearings  of  constant  surface  depth.  These 
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bearings  also  correspond  to  the  maximum  cross-range  slopes  of  the  surface,  where  we 
would  expect  to  see  the  strongest  amount  of  azimuthal  coupling. 
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Figure  3 1 .  Model  solutions  for  the  constant  slope  surface  perturbations  along  the  0  radial; 

(left)  using  the  basic  2D  interface  perturbations  model;  (right)  using  the  3D 

model. 
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Figure  32.  Model  solutions  for  the  constant  slope  surface  perturbations  along  the  180  radial; 

(left)  using  the  basic  2D  interface  perturbations  model;  (right)  using  the  3D 

model. 
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Figure  33.  The  wave  pattern  results  for  a  50  m  deep  source:  (left)  using  the  basic  2D 
interface  perturbations  model;  (right)  using  the  3D  model. 
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Figure  34.  Model  solutions  for  the  constant  slope  surface  perturbations  effect  on  a  50  m  deep 
source  waves:  (left)  using  the  basic  2D  interface  perturbations  model;  (right) 

using  the  3D  model. 
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Figure  35.  Difference  in  transmission  loss  values  at  a  constant  depth. 

2,  Downward  Sinusoidal  Surface  Perturbations 

The  second  environment  we  will  consider  for  comparison  between  the  fully  3-D 
results  versus  the  Nx2-D  approach  corresponds  to  the  single,  downward  sinusoidal 
surface  perturbation.  This  perturbation  produced  the  strongest  scattering  in  the  near-field, 
and  would  be  expected  to  produce  strong  differences  between  3-D  and  Nx-2D  in  the 
near-field,  as  well. 

Figures  36  through  39  present  the  results  from  the  two  approaches  for  this 
environment.  It  is  worth  noting  that  the  single  radial  results  presented  in  Figure  36  show 
some  marked  differences  in  the  field  computed  in  the  bottom  when  the  full  3-D  model  is 
employed.  This  is  due  to  high  vertical  angle  scattering  produced  by  the  computational 
algorithm.  The  transmission  loss  levels  of  these  differences  are  very  low,  however,  nearly 
120  dB  below  the  primary  energy  observed  in  the  water  column,  and  are  not  of  any 
physical  significance. 
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Figure  36.  Model  solutions  for  the  downward  sinusoidal  surface  perturbations  along  the  0 
radial:  (left)  using  the  basic  2D  interface  perturbations  model;  (right)  using  the  3D 

model. 
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Figure  37.  Model  solutions  for  the  downward  sinusoidal  surface  perturbations  along  the  180 
radial:  (left)  using  the  basic  2D  interface  perturbations  model;  (right)  using  the  3D 

model. 


As  before,  the  polar  plots  of  the  transmission  loss  field  at  a  constant  depth  of  30  m 
displayed  in  Figure  39  shows  little  apparent  difference  in  the  solutions.  In  an  effort  to 
highlight  the  differences,  Figure  40  shows  the  difference  in  transmission  loss  at  a 
constant  depth  in  polar  coordinates.  As  before,  the  differences  are  observed  to  expand  in 
range,  and  are  strongest  in  regions  where  the  cross-range  slope  of  the  surface  perturbation 
is  largest. 
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Figure  38.  The  wave  pattern  results  for  a  30  m  deep  souree;  (left)  using  the  basie  2D 
interface  perturbations  model;  (right)  using  the  3D  model. 
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Figure  39.  Model  solutions  for  the  downward  sinusoidal  surface  perturbations  on  a  30  m 
depth  source  waves:  (left)  using  the  basic  2D  interface  perturbations  model; 

(right)  using  the  3D  model. 
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Figure  40.  Difference  in  transmission  loss  values  at  a  constant  depth. 
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V.  CONCLUSIONS  AND  FUTURE  WORK 


Rough  surface  scattering  is  an  important  mechanism  for  introducing  signal 
variability  into  underwater  acoustic  systems  such  as  acoustic  communications.  The 
significance  of  three-dimensional  scattering  effects  from  two-dimensional  rough  surfaces 
has  been  investigated  by  other  researchers.  However,  whether  the  full  3-D  calculations 
constitute  a  significant  improvement  to  the  simpler,  Nx2-D  calculations  remains  unclear. 
The  goal  of  this  thesis  was  then  to  upgrade  an  existing  2-D  parabolic  equation  model  to 
incorporate  full  3-D  surface  scattering.  This  will  allow  future  studies  that  can  more 
directly  compute  the  specific  influence  of  3-D,  out-of-plane  scattering  from  2-D  rough 
surfaces. 

There  were  two  primary  advances  made  in  this  work  towards  the  stated  goal. 
First,  a  hybrid  implementation  for  the  parabolic  equation  model  was  formulated  and 
incorporated  into  a  version  of  the  Monterey-Miami  Parabolic  Equation  model.  This 
hybrid  method  consisted  of  the  standard  split-step  Fourier  algorithm  traditionally 
employed  in  MMPE,  and  also  a  finite  difference  method  utilizing  a  Pade  approximation 
for  the  cross-range  operator  associated  with  the  cross-range  scattering  of  the  rough 
surface.  The  solution  to  the  finite  difference  system  of  equations  was  accomplished 
through  a  standard  tri-diagonal  solver. 

The  second  advance  related  to  the  wavenumber  propagator  functions  employed  by 
the  3-D  model  in  cylindrical  coordinates.  Previous  work  assumed  the  azimuthal  coupling 
was  small  enough  that  a  standard  parabolic  equation  approximation  was  sufficient. 
However,  with  some  of  the  larger  surface  slopes  investigated  in  this  thesis,  it  was 
determined  that  a  wide  angle  parabolic  equation  approximation  was  necessary  in  order  to 
properly  account  for  the  cross-range  scattering.  Once  this  upgraded  propagator  was 
incorporated,  the  model  results  became  physically  realistic  for  all  of  the  rough  surfaces 
investigated. 

Five  different  “rough”  surfaces  were  defined  for  model  testing.  These  were  all 
simple,  non-fiat  boundaries  that  increased  in  complexity.  The  first  example  studied  was 
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simply  a  smooth,  sloped  surface.  The  resulting  field  appeared  consistent  with 
expectations.  In  order  to  validate  these  results,  the  model  was  recast  to  compute  the 
solution  for  a  flat  surface  with  a  sloping  bottom.  Through  an  environmental  rotation, 
these  results  should  be  identical.  A  comparison  of  the  model  results  confirmed  that  the 
sloping  surface  results  produced  valid  solutions. 

The  second  surface  was  defined  by  concentric,  sinusoidal  surface  displacements. 
Being  independent  of  bearing  angle,  this  surface  was  designed  to  test  the  model 
calculations  and  ensure  that  no  erroneous  out-of-plane  scattering  was  occurring  when 
none  was  expected.  The  model  results  were  consistent  with  this  expectation. 

The  third,  fourth,  and  fifth  surface  fluctuations  were  defined  in  terms  of 
horizontally  planar  sinusoidal  surface  perturbations  of  varying  types.  While  there  was  no 
formal  solution  to  compare  with  the  model  results,  the  general  behavior  was  consistent 
with  expectations,  where  scattering  was  observed  parallel  to  the  surface  wavefront. 

Finally,  model  results  for  two  of  the  surface  perturbations  were  computed  using 
both  the  3-D  version  developed  here,  as  well  as  the  Nx2-D  version  previously  developed. 
A  comparison  of  the  results  showed  the  dominant  differences  along  the  directions  with 
the  greatest  cross-range  surface  slopes,  as  expected. 

Future  work  could  entail  longer  range  solutions  in  order  to  search  for  more 
obvious  cross-range  coupling.  Using  sources  defined  by  individual  waveguide  modes 
could  also  be  employed  to  more  clearly  delineate  horizontal  refraction.  In  addition,  the 
model  results  examined  in  this  thesis  were  computed  using  only  a  single  frequency. 
Future  work  could  examine  broadband  signals,  and  thereby  compute  pulse  propagation 
travel  times,  which  could  also  provide  specific  evidence  of  horizontal  refraction  and  out- 
of-plane  scattering. 
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